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Abstract 

In this paper, we propose numerical methods for computing the boundary local 
time of reflecting Brownian motion (RBM) in M and its use in the probabilistic 
representation of the solution of the Laplace equation with the Neumann boundary 
condition. Approximations of the RBM based on a walk-on-spheres (WOS) and ran¬ 
dom walk on lattices are discussed and tested for sampling the RBM paths and their 
applicability in finding accurate approximation of the local time and discretization 
of the probabilistic formula. Numerical tests for several types of domains (cube, 
sphere, and ellipsoid) have shown the convergence of the numerical methods as the 
length of the RBM path and number of paths sampled increase. 

Key words: Reflecting Brownian Motion, Brownian motion, boundary local time, 
Skorohod problem, WOS, random walk, Laplace equation 


1 Introduction 


Traditionally numerical solutions of boundary value problems for partial differ¬ 
ential equations (PDEs) are obtained by using finite difference, finite element 
or boundary element methods with both space and/or time discretizations. 
This usually requires spatial mesh hue enough to ensure accuracy, which re¬ 
sults in considerable storage space requirement and computation time. More¬ 
over, the solution process is global, namely, the solutions of the PDEs have 
to be found together at all mesh points. However, in many scientihc and en¬ 
gineering applications, local solutions are sometimes all we need, such as the 
local electrostatic potential on a molecular surface where molecular binding 
activities are most likely to occur or the stress held at specihc locations where 
the materials are susceptive to failures. Therefore, it is of practical importance 
to have a numerical approach which can give a local solution of the PDEs at 
some locality of our choice. In the case of elliptic PDEs, this kind of local 
numerical method can be constructed using the well-known probabilistic rep¬ 
resentation and the associated Feynman-Kac formula raia. which relate ltd 
diffusion paths to the solution of the elliptic PDEs. By sampling the diffusion 
paths, the evaluation of the solution at any point in the domain can be done 
through an averaging process of the boundary (Dirichlet or Neumann) data 
under some given measure on the boundary. Moreover, this method avoids 
the expensive mesh generations required by mesh-based methods mentioned 
above na. 

Our previous work [1] , using the Feynman-Kac formula for the Laplace equa- 




tion with Dirichlet data, has produced a local method to compute the DtN 
(Dirichlet-to-Neumann) mapping for the Laplace operator. In this paper, we 
will focus on solving the following Neumann boundary value problem of the 
elliptic PDE using a probabilistic approach: 


{ An = /, on D 

= 0, ondD' 
on 


( 1 ) 


where D is a bounded domain in R^, A is the Laplace operator, / is a mea¬ 
surable function and 0 is a bounded measurable function on the boundary 

dD satisfying / (pda = / fdx. Equation (1) becomes the Laplace equation 
JdD Jd 

when / = 0, which is the subject of our work. 


The PDE ([T]) originates from either the Poisson equation for electrostatic 
potentials [19], an implicit time discretization of the heat equation or the 
momentum equation of the Naiver-Stokes equation with an additional lower 
term in the latter cases. Historically, Brownian motion (BM) has been used 
in solving PDEs due to its effectiveness and easy implementation regardless 
of dimensions jlO]. The well-known probabilistic representation can solve the 
elliptic equation with the Dirichlet boundary condition by using the first exit 
time td of BM, i.e.. 


u{x) = E%cP{xr^)) + 


r f{Xt)dt 

Jo 


( 2 ) 


In the above formula, only the values at the hitting positions on the boundary 
are used in the computation of the mathematical expectation (average) to 
obtain u{x). Taking the idea of killed Brownian motion [1] in junction with 
Monte Carlo methods, we can easily obtain an estimate of u{x). 


However, for the Neumann problem to be studied here, in contrast to the 
Brownian motion in (Ej), reflecting Brownian motion (RBM) will be needed 
to produce a similar probabilistic solution to ([T]). This theory has been devel¬ 
oped in [1] by employing the concept of the boundary local time whose one 
dimensional predecessor was introduced by Levy in [3|. In [1], the boundary 
local time of a one dimensional BM was extended to high dimensions and an 
explicit form, shown in ([2]), was obtained for domains with smooth boundaries. 
It should be noted that the boundary local time is related to the Skorohod 
equation [H] and plays a signihcant role in the theoretical development of the 
probabilistic approach to the Neumann problem. 


One-dimensional local time of Brownian motion has been studied by many 
authors rainiiEiiiii. For higher dimensions, similar results have been found 
by Brosamler [2]. Morillon [8] gave a modified Feynman-Kac formula for the 
Poisson problem with various boundary conditions, algorithms based on ran- 




dom walk on a grid have been proposed. However, numerical algorithms for 
computing local time in based on a rigorous probabilistic theory has not 
been done in the literature. It is the objective of this paper to obtain practical 
numerical methods for computing the local time of RBM in three dimensions 
and apply the resulting numerical methods to implement computationally the 
probabilistic representation for the Neumann problem. 

The rest of the paper is organized as follows. In section 2, we give some back¬ 
ground information on the Skorohod problem which is the key to the Neumann 
problem. In section 3, an explicit probabilistic solution to the Neumann prob¬ 
lem will be given. In section 4, a walk-on-spheres (WOS) method is reviewed 
and discussed for its application for RBM. In section 5, a numerical method, 
the WOS combined with a Monte Carlo method, is proposed for an approx¬ 
imation to the Neumann problem. Numerical results for cubic, spherical and 
ellipsoidal domains will be given in Section 6. Finally, we draw conclusions 
from our Monte Carlo simulations and discuss possible further work. 


2 Skorohod Problem, RBM and Boundary Local Time 


For the sake of completeness, we first give the dehnitions of Brownian motion 
and reflecting Brownian motion in R!^. 

Definition 1 Brownian motion: A Brownian motion B(t) = {Bi(t), B 2 {t),B^ 
in R'^ is a set of d independent stochastic processes with the following proper¬ 
ties: for I < i < d, 

(1) (Normal increments) Bi(t) — R(s) has a normal distribution with mean 
0 and variance t — s. 

(2) (Independence of increments) Bi(t) — Bi{s) is independent of the past, 
i.e., of Bu, 0 <u < s. 

(3) (Continuity of paths) Bi(t),t > 0 is a continuous function oft. 


Definition 2 Skorohod equation: Assume D is a bounded domain in R!^ with 
a boundary. Let f{t) be a (continuous) path in R^ with /(O) G D. A 
pair {(,t,Lt) is a solution to the Skorohod equation S{f]D) if the following 
conditions are satisfied: 

(1) f is a path in D; 

(2) Lit) is a nondecreasinq function which increases only when f G dD, 
namely, 

L{t)= f IaD{m)L{dsy, 

JO 


( 3 ) 



(3) The Skorohod equation holds: 


S(f\ D) : m = m -\f n(i(s))L(ds), (4) 

where n{x) stands for the outward unit normal vector at x & dD. 

Remark 3 In Definition 2, the smoothness constraint on D can be relaxed to 
hounded domains with boundaries, which however will only guarantee the 
existence o/(jl]). But for a domain D with a boundary, the solution will be 
unique. Obviously, {ft, Lt) is continuous in the sense that each component is 
continuous. 

If f{t) is replaced by the standard Brownian motion (BM) Bt, the correspond¬ 
ing ft will be a standard reflecting Brownian motion (RBM) Xt. Just as the 
name suggests, a reflecting BM (RBM) behaves like a BM as long as its path 
remains inside the domain D, but it will be reflected back inwardly along the 
normal direction of the boundary when the path attempts to pass through the 
boundary. The fact that Xt is a diffusion process with the Neumann boundary 
condition can be proven by using a martingale formulation and showing that 
Xt is the solution to the corresponding martingale problem with the Neumann 
boundary condition [1]. The result gives an intuitive and direct way to con¬ 
struct RBM from BM. This construction will be discussed in detail in Section 
5. 


Next we will review the concept of boundary local time L{t) for a RBM, which 
in a sense is a measure of the amount of time a RBM spends near the boundary 
and at the same time the frequency that a RBM hits the boundary. We have 
the following properties of L{t): 


(a) R is the unique continuous nondecreasing process that appears in the 
Skorohod equation (|1]) [1]; 

(b) It measures the amount of time the standard reflecting Brownian motion 
Xt spending in a vanishing neighborhood of the boundary within the 
period [0,t]. If D has a boundary, then 


L{t) = lim 


fpD,(X,)ds 

6 


(5) 


where D^ is a strip region of width e containing dD and D^ G D. This 
limit exists both in and P^-a.s. for any x E D] 

(c) L{t) is a continuous additive functional (CAF) |2] which satishes the 
additivity property n HUS] HU: At+s = As + At{6s). Here Og denotes the 
shift operator along the paths. The additivity property of L{t) can be 
seen as follows: 



From the definition in ([5]), we have 


L{t + s) 
L{s) 
L{t) o Os 


Io^{X,)dT 

fo^lDAk)dT 


lim 

e->0 

lim 

e-i^O 

AlnAXr+s)dT ^ r,^^InAXr)dr 

e —>-0 f e —>-0 e 


therefore, 


F(t + s) — -h(t) + L(t) o 9g^ 

which shows that L(t) satishes the additivity property. 


For one-dimensional case, much existing literature devoted to the study of 
local times for Brownian motion and more general diffusion processes. It is well 
known that for one dimensional Brownian motion starting from the origin, the 
local time L{t) of RBM and maxi?(s) have the same distribution as stochastic 

s<t 

processes. Hence, valuable properties of RBM can be drawn by just observing 
maxR(s). But this is not true in general in higher dimensions. However, we 

s<t 

have the following explicit formula for L{t) derived in [T], 

m = ^f^l9D{Xg)^S, ( 6 ) 

where the the right-hand side of ([6]) is understood as the limit of 


n-l 

i=\ 

where A = {Aj} is a partition of the interval [0,f] and each Aj is an element 
in A. We will discuss the implementation of both (|5]) and (jb]) in Section 5. 



max/a£)(A^ 

seAi 


Aj|, max|Aj| —)■ 0, 


(7) 


3 Neumann Problem 


We will consider the elliptic PDF in with a Neumann boundary condition 









When the bottom of the spectrum of the operator A/2 + g is negative a 
probablistic solution of (E]) is given by 


u{x) = -E^ / eq{t)(j){Xt)L{dt) 
2 Uo 


(9) 


where Xt is a RBM starting at x and eq{t) is the Feynman-Kac functional |T] 

eq{t) = exp [ q{Xs) ds 
Jo 


From the dehnition of the local time in (|5]), we have the following approxima¬ 
tion for small e 

/Pa (A',)* 


L(t) 

Plugging ffTOj) into ([9]), we have 

, . 1 


u[x) ^ 


POO pt-\-dt 

e,(t)0(W) lDXXs)ds 


( 10 ) 


( 11 ) 


The solution defined in ([9]) should be understood as a weak solution for the 
classical PDF (E]). The proof of the equivalence of (Ej) with a classical solution 
is done by using a martingale formulation |T]. If the weak solution satishes 
some smoothness condition nH. it can be shown that it is also a classical 
solution to the Neumann problem. This formula is the basis for our numerical 
approximations to the Neumann problem ([H]). To compute the expectation in 
the formula, we rely on Monte Carlo random samplings to simulate Brownian 
paths and then take the average. 


In the present work, as we only consider the Laplace equation where g = 0, 
therefore. 


u{x) 


—E^ 

2e 


ps-\-dt 


0(W) / lDXXs)ds , (12) 

/O Js 

and we will show how this formula is implemented with the Monte Carlo and 
WOS methods in section 5. 


Remark 4 Comparing with formula (E]), we find that the probabilistic solu¬ 
tions to the Laplace operator with the Dirichlet boundary condition has a very 
similar form (referring to flTT)) ). In the Dirichlet case, killed Brownian paths 
were sampled by running random walks until the latter are absorbed on the 
boundary and u{x) is evaluated as an average of the Dirichlet values at the 
first hitting positions on the boundary, namely, u{x) = E^ [fi^Xrj^)] where cf is 
the Dirichlet boundary data. On the other hand, for the Neumann condition, 
while u{x) is also given as a weighted average of the Neumann data at hitting 











positions of RBM on the boundary, the weight is related to the boundary local 
time of RBM. This is a noteworthy point when we compare the probabilistic so¬ 
lutions of the two boundary value problems and try to understand the formula 
in dH]). 


4 Method of Walk on Spheres (WOS) 


Random walk on spheres (WOS) method was hrst proposed by Muller [7], 
which can solve the Dirichlet problem for the Laplace operator efficiently. 
Here we will hrst briehy review this method and then show how it can be 
adapted for RBM and the Neumann problem. 


For a general linear elliptic problem with a Dirichlet boundary condition, 


n O n 02 

L{u) = ^k{x)-— + ^ = /(^)>^ e D, 

OX, UXiUXj 

u\dD = 4>{x),x e dD. 

The probabilistic representation of the solution is (116] [H]) 


U{x) = E^{(f){Xro)) + 


f{Xt)dt 


where Xt{w) is an Ito diffusion defined by 

dXt = b{Xt)dt + a{Xt)dBt, 


(13) 


(14) 


(15) 


and Bt{w) is the Brownian motion, [ojj] = -a(x)a^(x). 


The expectation in ffTT|) is taken over all sample paths starting from x and 
td is the first exit time for the domain D. This representation holds true for 
general linear elliptic PDFs. For the Neumann boundary condition, similar 
formulas can be obtained [S]. However different measures on the boundary 
dD will be used in the mathematical expectation. 


In order to illustrate the WOS method for the Dirichlet problem, let us con¬ 
sider the Laplace equation where / = 0, = 6ij and 5, = 0 in flT3l) and the 

Ito diffusion is then simply the standard Brownian motion with no drift. The 
solution to the Laplace equation can be rewritten in terms of a measure /if, 
defined on the boundary dD, 


'll(x) = E‘(tl>(Xr„)) = 


( 16 ) 





where /i|) is the so-called harmonic measure defined by 


/i^(F) = P" e P}, P c ap,X e P. (17) 


It can be shown that the harmonic measure is related to the Green’s function 
for the domain with a homogeneous boundary condition, 


-^g{.x,y) = 5{x - y), x e D, 
g{x, y) = 0, X E dD 


(18) 


By the third Green’s identity, 


u{x) 



u{y) 


dg{y,x) 

dn 


9{y,x)^{y) 


dSy, 


(19) 


and using the zero boundary condition of we have 

u{x) = f u{y)^^^^^dSy. ( 20 ) 

JdD on 

Thus, the hitting probability y]j{[y,y + dSy]) is equivalent to p{x, y)dSy. Gom- 
paring flTOj) with fl20|) . we can see that 


P(x,y) = 


dn,. 


( 21 ) 


For instance, the Green’s function for a ball for this purpose is given as 

1 1 


g{x,y) = 


+ 


4:7i\x — y\ 47r|x —?/*|’ 
where y* is the inversion point of y with respect to the sphere [1]. 


( 22 ) 


If the starting point x of a Brownian motion is at the center of a ball, the 
probability of the BM exiting a portion of the boundary of the ball will be 
proportional to the portion’s area. It is known that all sample functions of 
Brownian motion processes starting in the domain intersects the boundary 
dD almost surely [7]. Therefore, sampling a Brownian path by drawing balls 
within the domain, regardless of how the path navigates in the interior before 
hitting the boundary, can significantly reduce the path sampling time. To 
be more specific, given a starting point x inside the domain P, we simply 
draw a ball of largest possible radius fully contained in P and then the next 
location of the Brownian path on the surface of the ball can be sampled, using 
a uniform distribution on the sphere, say at xi. Treat xi as the new starting 
point, draw a second ball fully contained in P, make a jump from xi to X 2 on 
the surface of the second ball as before. Repeat this procedure until the path 
hits a absorption e-shell of the domain [5]. When this happens, we assume 
that the path has hit the boundary dD (see Figure 1(a) for an illustration). 











(b) WOS (with a maximal 
step size for each jump) within 
the domain 


Fig. 1. Walk on Spheres method 



Fig. 2. A e-region for a bounded domain in 
Next we define an estimator of ffT^ by 

1 ^ 

~ ( 23 ) 

i=i 

where N is the nnmber of Brownian paths sampled and Xi is the hrst hitting 
point of each path on the bonndary. Using a jnmp size (radius of the ball) 5 
on each step for the WOS, we expect to take 0(1/5^) steps for a Brownian 
path to reach the boundary [6j. To speed up, maximum possible size for each 
step would allow faster first hitting on the boundary. Most of the numerical 
results in this paper will use the WOS approach as illustrated in Figure 1(b). 


5 Numerical Methods 


5.1 Simulation of reflecting Brownian paths 


A standard reflecting Brownian motion path can be constructed by reflecting 
a standard Brownian motion path back into the domain whenever it crosses 
the boundary. So in principle, the simulation of RBM is reduced to that of 
BM. 

It is known that standard Brownian motion can also be constructed as the 
scaling limit of a random walk on a lattice so we can model BM by a random 









Fig. 3. WOS in the e-region. BM path hits xi in e-region for the first time. Then the 
radius of sphere is changed to Ax, the path continues until it arrives at X 2 whose 
distance to dD is smaller or equal to Ax. Then the radius of the ball is enlarged to 
2Ax so that the path has a chance to run out of the domain at X3. If that happens, 
we pull back X 3 to X 4 which is the closest point to X 3 on the boundary. Record 
4>{xi), and continue WOS-sampling the path starting at X4. 

walk with proper scale (see Appendix for details). However, it turns out that 
the WOS method is the preferred method to simulate BM for our purpose |9] 
(see Remark [6] for details). As mentioned before, a e-shell is chosen around 
the boundary as the termination region in the Dirichlet case. Here we follow 
a similar strategy by setting up a e-region but allowing the process Xt to 
continue moving after the latter reaches the e-region instead of being absorbed. 

Figure 2 shows a strip region with width e near the boundary is identihed for 
a bounded domain. In a spherical domain, the e-region is simply an e-shell 
near the boundary of width e. Denote M^{D) as the e-region and I{D) as the 
remaining interior region D\M^{D). 

Recall the discussion of the WOS in the previous section. For a BM starting 
at a point x in the domain, we draw a ball centered at x, the Brownian path 
will hit the spherical surface with a uniform probability as long as the ball 
does not overlap the domain boundary dD. The balls are constructed so that 
the jumps are as large as possible by taking the radius of the ball to be the 
distance to the boundary dD. We repeat this procedure until the path reaches 
the region M^{D). Here, we continue the WOS in M^{D) but with a hxed 
radius Ax much smaller than e. In order to simulate the path of RBM, at 
some points of the time the BM path may run out of the domain. For this 
to happen, the radius of WOS is increased to 2Ax when the path is close 
to boundary at a distance less than Ax. In this way, the BM path will have 
a chance to get out of the domain, and when that happens, we then pull it 
back to the nearest point on the boundary along the normal of the boundary. 
Afterwards, the BM path will continue as before. 




Fig. 4. A RBM path with a cube in 

In summary, a reflecting Brownian motion path is simulated by the WOS 
method inside D. Once it enters the e-region M^{D), the radius of WOS 
changes to a fixed value, either Ax or 2Ax, depending on its current distance 
of the Brownian particle to the boundary. Once the path reaches a point on 
the boundary after the reflection, the radius of WOS changes back to Ax. 
Figure 3 illustrates the movement of RBM in the e-region M^{D). As time 
progresses, we expect the path hits the boundary at some time instances and 
lies in either I{D) or M^{D) at others. A RBM path is shown in Figure 4 
within a cube of size 2. 


5.2 Computing the boundary local time L{t) 


Two equivalent forms of the local time have been given in ([5]) and dH]). Here we 
will show how the e-region for the construction of the RBM in Fig. 3 can also 
be used for the calculation of the local time. When the e-region is thin enough, 
i.e. e <C 1, an approximation of ([5]) is given in flTOl) . which is the occupation 
time that RBM Xg sojourns within the e-region during the time interval [0, t]. 
A close look at flTOl) reveals that only the time spent near the boundary is 
involved and the specific moment when the path enters the e-region has no 
effect on the calculation of L{t). 

Suppose X E D is the starting point of a Brownian path, which is simulated 
by the WOS method. Once the path enters the e-region, the radius of WOS 
is changed to Ax or 2Ax. It is known that the elapsed time At for a step 
of a random walk on average is proportional to the square of the step size, 
in fact. At = {Ax)‘^/d,d = 3 when Ax is small (see Appendix), which also 
applies to WOS moves (See Remark 6 for details). Therefore, we can obtain 
an approximation of the local time L{dt) by counting the number of steps the 




Fig. 5. Boundary local time (1241) increases when the path runs into the region M^{D). 
The insert shows the piecewise linear profile of the local time path with flat level 
regions. The path of L{t) is a nondecreasing function and increases at a rate lower 
than 0{NT). NT is the length of the path. 

path spent inside M^{D) mnltiplied by the time elapsed for each step, i.e. 


L{dt) = L{tj — tj-i) 


r^j 

Jn 


lDSXs)ds 


= {nt. - 


(Ax)" 

3e 


(24) 


where n*. —nt _^ is the nnmber of steps that WOS steps remain in the e-region 
dnring the time interval [tj-i, tj]. Fignre 5 gives a sample path of the simulated 
local time associated with the RBM in Figure 4. 


Remark 5 (Alternative way to compute local time L{t) ) From ([6]), the local 
time increases if and only if the RBM path hits the boundary, which implies 
that the time before the path hits the boundary makes no contribution to the 
increment of the local time. Thus, a WOS method with a changing radius can 
also be used with ([6]). Specifically, we divide the time interval [0,f] into to 
N small subintervals of egual length. In each the Brownian path will 

move 2Ax or Ax with the WOS method when the current path lies within a 
distance less or more than Ax to the boundary. If the path hits or crosses the 


boundary within then L{t) will increase by \J'n 


i+l 


Remark 6 (Approximating RBM by WOS or random walks on a lattice - 
a comparison) There are two ways to find approximation to Brownian paths 
inside the region M^{D) and construct their reflections once they get out of the 
boundary. One way is by using the WOS approach as described in Section 5.1. 
The other is in fact to use a random walk on a lattice inside MflD). In the 
second approach, as illustrated in Fig. 6, a grid mesh is set up over MflD) and 
the random walk takes a one-step walk on the lattice until the path goes out of 
the domain and then it will be pushed back to the nearest lattice point inside 
MflD). And the elapsed time for a Ax walk is on average {Axfl/3 as shown in 














Fig. 6 . Random walks on the e-region. A BM path hits xi E M^{D) by the WOS 
method. Replace xi by the nearest grid point x'l. Then several steps of random 
walks will make a path as X 2 —>■ ica —)• X 4 . Since X 4 ^ T), we push it back along 
the normal line (dash arrow) to X 4 then replace it by the closest grid point within 
domain (solid arrow) X5. Here path crosses the boundary at X4 E dD. Then continue 
the random walk as usual at xq. 

the Appendix. The boundary local time L(t) can be still calculated as in (Elj). 
The problem with this approach is that a Brownian motion actually should 
have equal probability to go in all directions in the space while a random walk 
on the lattice only considers six directions in . This limitation was found 
in our numerical tests to lead to insufficient accuracy in simulating reflecting 
Brownian motions for our purpose. 

Meanwhile, the WOS method in the e-region MflD) has a fixed radius Ax, 
which enables us to calculate the boundary local time by since the elapsed 
time of a Ax move in on average still remains to be (Ax)^/3. This conclu¬ 
sion can be heuristically justified by considering points on the sphere are the 
linear combination of the directions along the three axes, which implies that 
the average time that the path hits the sphere with radius Ax should also be 
the same. As discussed before, if the path comes within a distance very close 
to the boundary, say less than Ax, the radius of the WOS method is increased 
to 2Ax so that it will have a chance to run out of the domain and then be 
pushed back to the nearest point on the boundary to affect a hit of the RBM 
on the boundary. 


5.3 Probabilistic representation for the Neumann problem 


Finally, with the boundary local time of RBM available, we can come to the 
approximation of the Neumann problem solution u{x) using the probabilistic 
approach (IT^ . First of all, we will need to truncate the inhnite time dura- 

















tion required for the RBM path Xt in ffT^ to a hnite extent for computer 
simulations. The exact length of truncation will have to be numerically deter¬ 
mined by increasing the length until a convergence is conhrmed (namely, the 
approximation to u{x) does not improve within a prescribed error tolerance 
between two different choices of truncation times under same number of sam¬ 
pled paths). Assume that the time period is limited to from 0 to T, then by a 
Monte Carlo sampling of the RBM paths, an approximation of flT^ will be 


u{x) 


1 

Ye 


N V 
i=i Ro 



lD.(Xi)da 


(25) 


where XI,i = are stochastic processes sampled according to the law 

of RBM. 

Next, let us see how the RBM can be incorporated into the representation 
formula once its path is obtained. 

Associate the time interval [0,T] with the number of steps NT of a sampling 
path, NT will give the total length of each path. Then, the integral inside the 
square bracket in fl2^ can be transformed into j 

NT / t. \ 

E ’t’N )-f®(A'),) / /d. (A))* , (26) 

j>=l \ “'b-1 / 

where j' stands for the j'—th step the WOS method has taken, and j indicates 
a step for which Xl^ G dD. 

As the integral in fl2BD is in fact the occupation time as shown in fl^TD . 
becomes 

E (^(Xl)IMX;,)(n,, - ■ (27) 


As a result, an approximation to the PDE solution u{x) becomes 


1 ^ 




NT 




(Ax) 


f=i 


(28) 


Theoretically speaking, e should be chosen much larger than Ax. Here, we 
take e = kAx, fc > 1 is an integer, which will increase as Ax vanishes to zero. 








Then, fl28|) reduces to 


uix = 


1 ^ 

^—y 

2kAx ^ 


NT 


Eh'R)WAl)K 


(Ao:) 


i=i 


Ax ^ 


NT 


i=i 


(29) 


which is the hnal numerical algorithm for the Neumann problem. In the fol¬ 
lowing we present the general implementation of this numerical algorithm. 


Let X be any interior point in D where the solution u{x) for the Neumann 
problem is sought. First, we define the e-region M^{D) near the boundary. For 
each one of N RBM paths, the following procedure will be executed until the 
length of the path reaches a prescribed length given by NT ■ Ax: 


(1) If X ^ M^{D), predict next point of the path by the WOS with a maximum 
possible radius until the path locates near the boundary within a certain 
given distance e, say e = 5Ax (hit the e-region M^{D)). If x G M^{D), 
l{ti) = 1; otherwise, l{ti) = 0. Here l{t) is the unit increment of L{t) at 
time t. 

(2) If X G use the WOS method with a hxed radius Ax to predict the 

next location for Brownian path. Then, execute one of the two options: 


Option 1. If the path happens to hit the domain boundary dD at x*., 
record 0(xtj. 


Option 2. If the path passes crosses the domain boundary dD^ then pull 
the path back along the normal to the nearest point on the boundary. Record 
the Neumann value at the boundary location. 


Due to the independence of the paths simulated with the Monte Carlo method, 
we can run a large number of paths simultaneously on a computer with many 
cores in a perfectly parallel manner, and then collect all the data at the end of 
the simulation to compute the average. Algorithm 1 gives a pseudo-code for 
the numerical realization of implementing the WOS in both I{D) and M^{D) 
regions. 








Data: Select integers N and NT, a starting point Xq G D, step size h and 
e-region M^{D) near the boundary. 

Output An approximation of u{Xo). 

Initialization L[NT],v[NT],u[N], X = Xq, i ■(— 1 and j ■(— 1; 

While i < N do 
Set Si = 0. 

While j < NT do 

If X G I{D) then /* If the path has not touched the e-region */ 
Set L[j] ■(— 0; /*Increment of local time at each step. */ 

Set r d{X, dD)-, /* Find the distance to the boundary */ 

Randomly choose a point Xi on B{X,r) then set X ■(— Xi. 

Else /* The path enters the e-region */ 

L[j] 1; /*local time increases */ 

Set r G- h (2h); /* If d{X, dD) > h or =0 (0 < d{X, dD) < h) */ 

Randomly choose a point Xi on B(X,r) then set X ■(— Xi. 

If X ^ D, then 

Find Xj to be the nearest point on dD to X and pull X back 
onto dD at Xj] 

Set X ^Xj] 

Set v[j] ^ (j){Xj) 

End 

End 

3^3 +1 ; 

End 

count ■(— 0 ; 

For k=l:NT 
count count -|- L[k\, 

If v[k] ~= 0 then 
u[i\ u[i\ -|- 0(Xfc)-count; 
count ■(— 0 ; 

End 

i ^ i + l] 

End 

N 

Return u{Xq) = h'^u[k]/N/{6k) 

k=l 


Algorithm 1: The algorithm for the probabilistic solution of the Laplace 
equation with the Neumann boundary condition 





6 Numerical results 


In this section, we give the numerical results for the Neumann problem in 
cubic, spherical and ellipsoid domains. 

To monitor the accuracy of the numerical approximation of the solutions, we 
select a circle inside the domain, where the solution of the PDE u{x) will be 
found by the proposed numerical methods, dehned by 

{(x, y, zY = cos 9i sin 02, r sin 6i sin 02, r cos O 2 Y} (30) 

with r = 0.6, 0i = 0 : A; • 27r/30 : 27r, 62 = vr/4 with /c = 1,..., 15. In addition, 
a line segment will also be selected as the locations to monitor the numerical 
solution, the endpoints of the segment are (0.4, 0.4, 0.6)^ and (0.1, 0,0)^, re¬ 
spectively. Fifteen uniformly spaced points on the circle or the line are chosen 
as the locations for computing the numerical solutions. 

The true solution of the Neumann problem ([1]) with the corresponding Neu¬ 
mann boundary data is 


u{x) = sin3xsin4|/ -|- 5. 


(31) 


In the hgures of numerical results given below, the blue curves are the true 
solutions and the red-circle ones are the approximations. The numerical solu¬ 
tions are shifted by a constant so they agree with the exact solution at one 
point as the Neumann problem is only unique up to an arbitrary additive 
constant. “Err”indicates the relative error of the approximations. 


6 .1 Cube domain and test on the length of the path 


A cube domain of size 2 is selected to test the choice of the number of paths 
and the length of the paths (truncation time duration T) in the numerical 
formula fl 2 ^ . 

The step-size Ax = 0.0005 is used as the radius of the WOS inside the e-region 
M^{D), namely, the step-size of the random walk approximation of the RBM 
near the boundary. The number of paths is taken as N = 2e5. Two choices for 
the path length parameter NP = 2.7e4 and NP = 3e4 are compared to gauge 
the convergence of the numerical formula fl 2 ^ in terms of the path truncation. 
Figures 7 and 8 shows the solution and the relative errors in both cases, which 
indicates that NP = 3e4 will be sufficient to give an error around 5% as shown 
in Fig. 8 . 




(a) e = 6Ax, Err = 10.71% (b) e = 7Ax, Err = 12.19% 

Eig. 7. Cubic domain: number of paths N = 2e5, and number of steps for each path 
NP = 2.7e4. (Left) Solution on the circle defined in (1,11 p . (right) solution on a line 
segment. 




(a) e = 6Ax, Err = 5.35% (b) e = 7Ax, Err = 5.85% 

Fig. 8. Cubic domain: number of paths N = 2e5, and number of steps for each path 
NP = 3e4. (Left) Solution on the circle defined in ()3ip . (right) solution on a line 
segment. 

In the rest of the numerical tests, we will set the number of path N = 2e5, 
and number of steps for each path NP = 3e4. 


6.2 Spherical domain 


The unit ball is centered at the origin. We set Nx = 0.0005 and adjust e, 
similar numerical results are obtained as in the case of the cube domain. Here, 
the reflected points of Brownian path are the intersection of the normal and 
the domain. Though Figure 8(b) shows some oscillations in the middle, the 
overall approximation are within a relative error around 5.85%. 


6.3 Ellipsoid domain 


The ellipse with axis lengths (3, 2, 1) is centered at the origin. We set Ax = 
0.0004. The numerical results along the circle behave better than those along 








(a) e = 5Ax, Err = 5.13% (b) e = 5Ax, Err = 4.03% 

Fig. 9. Spherical domain: number of paths N = 2e5, and number of steps for each 
path NP = 3e4. (Left) Solution on the circle defined in (lOTT) . (right) solution on a 
line segment. 




(a) e = 5Ax, Err = 5.75% (b) e = 5Ax, Err = 5.12% 

Fig. 10. Ellipsoid domain: number of paths N = 2e5, and number of steps for each 
path NP = 3e4. Left: solution on the circle defined in (|3ip . Right: solution on a 
line segment. 

the line segment, especially along the tail section of the latter (Figure 10(b)), 
which he closer to the origin 0. 


7 Conclusions and discussion 


In this paper we have proposed numerical methods for computing the lo¬ 
cal time of reflecting Brownian motion and the probabilistic solution of the 
Laplace equation with the Neumann boundary condition. Without knowing 
the complete trajectories of RBM in space, we are able to use the WOS to sam¬ 
ple the RBM and calculate its local time, based on which a discrete probabilis¬ 
tic representation ( 12 ^ was obtained to produce satisfactory approximations 
to the solution of the Neumann problem at one single point. Numerical results 
validated the stability and accuracy of the proposed numerical methods. 

In addition, random walk on a lattice was also investigated as an alternative 
way to sample RBM. However, numerical experiments show that the numerical 






results are inferior to those obtained by the WOS method. A possible reason 
is that formula (|5]) for the local time is valid for a smooth path while a random 
walk approximation of the the Brownian path contains inherent errors. 

The local time can also be computed by a mathematically equivalent formula 
(|6]), for which the implementation is discussed briefly in section 5.2. Again the 
numerical results based on (I6]) are inferior to those obtained using the original 
limiting process of Levy in [3] . This fact we believe may result from the time 
discretization error of Brownian paths especially when long time truncation is 
employed in the probabilistic representation. 

Various issues affecting the accuracy of the proposed numerical methods re¬ 
main to be further investigated, such as the number of random walk or WOS 
steps and the truncation of duration time T for the paths, the choice of the 
thickness for the e-region, the size of Ax for the lattice, etc. In theory, the 
larger the truncation time T, the more accurate is the probabilistic formula 
for the Neumann solution. However, for a hxed spatial mesh size Ax, long time 
integration will result in the accumulation of time discretization error for the 
Brownian pathes, thus leading to the degeneracy of the numerical solutions as 
verihed by our numerical experiments. Meanwhile, Revesz [18] have proposed 
some approximations of local time by other stochastic processes in the case 
of a half line. We conjecture such results may still hold in higher dimensions 
and progress in this direction will shed light on how to improve the numerical 
procedures proposed in this paper. 


8 Appendix 


If the random walk on a lattice as in Fig. 11 
BM, a relationship between At and Ax in 
be 


is to converge 
will be needed 


to a continuous 
and is shown to 


At = 


(Ax)^ 

3 


(32) 


The following is a proof for this result (See [I2] for reference). The density 
function of standard BM satishes the following PDF [1] 


dp 

m 


1 

2 


l^xP(t,x,y) . 


(33) 


By using a central difference scheme and changing p to u, equation (j33il be- 



(x,y,z+l) 



Fig. 11. Central difference scheme in 


comes 


At 


{Ax) 


(34) 


Reorganizing and letting A = At/(2(Aa;)^) , we have 

= ^'^i+i,j,k+^'Vi-i,j,k+^'^i,j+i,k+^'^i,j-i,k+^'^i,j,k+i+^'^i,j,k-i+i^~^^)'^i,j,k ) 

^ (35) 

By setting A = -, we have 
6 

'^Sk = g^i+i,i,fc + ^'^i-i,j,k + ^'^i,j+i,k + -^'^i,j-i,k + '^'^i,j,k+i + ^'^i,j,k-i- (36) 


6 


6 


6 


6 


6 


For the initial condition 0, we have 


where 


<d= E 

\/=l 


(-h,0,0)^, 

prob 

(h,0,0)^ 

prob 

(0,h,0)^ 

prob 

1 

o 

prob 

o' 

o"' 

prob 

Sh 

1 

o' 

o'" 

prob 


1 

6 

1 

6 

1 

6 

1 ’ 

6 

1 

6 

1 

6 


(37) 


vi = s 


(38) 



and 


E^ = 


i=i 


^ — n + 2i' + i \ 

— n + 2j' + j h. 
^ — n + 2k' + kj 


(39) 


Let r]i = {xi,yi,zi)^, then 


xi = < 


-h, prob = - 
6 

h, prob = - , 
6 
2 

0, prob = - 


(40) 


for each 1. We known that yi, zi have the same distribution as xi. 


Notice that the covariance between any two of xi, yi, zi is zero, i.e. E{xiyi) = 0, 

n n n n 

E{yizi) = 0 and E{xiZi) = 0. So E{'^xi'^yi) = 0, E{'^yi'^zi) = 0 and 

i=l i=l i=l i=l 

n n 

EC^Xi'^zi) = 0. According to the central limit theorem, we have 

i=l i=l 


n / ^/j2\ 

y^ xi = N 1^0, n —)■ cxD. (41) 


n n 

The same assertion holds for '^yi and zi. 

i=\ i=\ 


Since A 

n 


At 

2{Axy 


-, then 
6 


= 3k and hence 


nh^ 


n 


n 


xz ~ A^(0, t) as n —)■ oo. So are yi and zi. 
2=1 2=1 2=1 


nk = t. Therefore 


n n n 

Recall that the covariance between any pair of '^xi, '^yi, and '^zi is zero, 

2=1 2 = 1 2 = 1 

n n n 

that and '^^zi are independent normal random variables. Hence, 

2 = 1 2=1 2 = 1 


C i'^jf ^ V Vl 


1=1 


^ — n + 2i' + i \ 

— n + 2j' + j h = 

— n + 2k' + kj 


M 

2 = 1 


Y.yi 

2 = 1 
n 

zi 

vs yj 


>4 


-H ^ -^olr 


(2ir«)S/2 


(42) 



and 


^ JJJr> {2Rtrn 


e ^ 2 r“— (p{j£)dx, (43) 


which coincides with the density function of the S-d standard BM. 

T 1 • 1 1 • A (Ax)^ r-^ dx ^ , 

in conclusion, when —-—— = i.e. At = - or vat = the central 

2 (Aa;)^ 6 3 V3 

difference scheme converges to the standard BM in 3-d. Generally, the result 
can be extended to d-dimensional Euclidean space and the result will be At = 
{Axf 

d ' 
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